# Load necessary libraries
library(AER)
library(ltm)
library(MASS)
library(broom)
library(dplyr)
library(psych)
library(Matching)
library(mediation)
library(dotwhisker)
library(modelsummary)
library(interactions)
source("process.R")

# Load data
data = readRDS("data.ii.rds")

# Cronbach Alpha results (see footnote 6 on page 12 )
print("Cooperative Cronbach Alpha")
cronbach.alpha(data.frame(data$coop.1, data$coop.2, data$coop.3, data$coop.4, data$coop.5))
print("Militant Cronbach Alpha")
cronbach.alpha(data.frame(data$mili.1, data$mili.2, data$mili.3, data$mili.4.R, data$mili.5))
print("Isolationist Cronbach Alpha")
cronbach.alpha(data.frame(data$isol.1, data$isol.2, data$isol.3.R, data$isol.4, data$isol.5))

# Descriptive Statistics (see Analysis Section)
print("Age Descriptive Statistics")
describe(data$age)
print("Male Participant Percentage")
prop.table(table(data$is.male))
print("Education Percentage")
prop.table(table(data$educ)) 
print("News Descriptive Statistics")
describe(data$news)
print("Ideology Descriptive Statistics")
describe(data$ideo) 
print("Cooperative Descriptive Statistics")
describe(data$coop)
print("Militant Descriptive Statistics")
describe(data$mili)
print("Isolationist Descriptive Statistics")
describe(data$isol)
print("Cyber Knowledge Descriptive Statistics")
describe(data$cybr.know)
print("Political Knowledge Descriptive Statistics")
describe(data$polt.know)
print("Cyber Response Outcome Descriptive Statistics")
describe(data$resp.cybr.s)
print("Conventional Response Outcome Descriptive Statistics")
describe(data$resp.conv.s)
print("Interest Mediator Descriptive Statistics")
describe(data$intr)
print("Identity Mediator Descriptive Statistics")
describe(data$idnt)
print("Reputation Mediator Descriptive Statistics")
describe(data$rept)
print("Morality Mediator Descriptive Statistics")
describe(data$morl)

print("T-Test Between Outcome Means")
t.test(data$resp.cybr.s, data$resp.conv.s, alternative = "two.sided") 

print("Correlation Analysis Between Outcome Means")
data %>% select(intr, idnt, rept, morl) %>% cor()

# Causal Analysis Defend
linr.defend.treat = lm(defend ~ commitment + casualties + dyad, data = data)
linr.defend.media = lm(defend ~ commitment + casualties + dyad + morl + intr + rept + idnt, data = data)
linr.defend.nmedia = lm(defend ~ commitment + casualties + dyad + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.defend.full = lm(defend ~ commitment + casualties + dyad + morl + intr + rept + idnt + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)

# See Online Appendix B

print("Defend - Treatment Model")
summary(linr.defend.treat)
print("Defend - Treatment + Mediator Model")
summary(linr.defend.media)
print("Defend - Full Model")
summary(linr.defend.full)

# Dot Whisker Plot (Figure 1)

dwplot(list(linr.defend.nmedia, linr.defend.full),
       vline = geom_vline(
         xintercept = 0,
         colour = "grey60",
         linetype = 2
       ),
#       dot_args = list(aes(shape = model)),
#       whisker_args = list(aes(linetype = model))
) %>%
  relabel_predictors(
    c(
      dyadYes = "Dyad",
      casualtiesYes = "Casualties",
      commitmentYes = "Commitment",
      morl = "Morality",
      intr = "Interest",
      idnt = "Identity",
      rept = "Reputation",
      ideo.c = "Ideology",
      news.c = "News",
      cybr.know.c = "Cyber Knowledge",
      polt.know.c = "Political Knowledge",
      coop.c = "Cooperative Internationalism",
      is.maleYes = "Male",
      age.c = "Age"
    )
  ) +
  xlab("Defend") + ylab("") +
  theme_bw() +
  scale_colour_grey(
    start = .2,
    end = .6,
    name = "Defend Models",
    labels = c('Without Mediators', 'With Mediators')
  ) +
  scale_fill_grey(
    start = .2,
    end = .6,
    name = "Defend Models",
    labels = c('Without Mediators', 'With Mediators')
  )

# Causal Analysis Mediators

linr.morl.mod = lm(morl ~ commitment * casualties + commitment * dyad + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.intr.mod = lm(intr ~ commitment * casualties + commitment * dyad + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.rept.mod = lm(rept ~ commitment * casualties + commitment * dyad + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.idnt.mod = lm(idnt ~ commitment * casualties + commitment * dyad + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)

# See Online Appendix E

print("Morality Model")
summary(linr.morl.mod)
print("Interest Model")
summary(linr.intr.mod)
print("Reputation Model")
summary(linr.rept.mod)
print("Identity Model")
summary(linr.idnt.mod)

# Dot Whisker Plot (Figure 2)

dwplot(list(linr.morl.mod, linr.intr.mod, linr.rept.mod, linr.idnt.mod),
       vline = geom_vline(
         xintercept = 0,
         colour = "grey60",
         linetype = 2
       )
) %>%
  relabel_predictors(
    c(
      dyadYes = "Dyad",
      casualtiesYes = "Casualties",
      commitmentYes = "Commitment",
      ideo.c = "Ideology",
      news.c = "News",
      cybr.know.c = "Cyber Knowledge",
      polt.know.c = "Political Knowledge",
      coop.c = "Cooperative Internationalism",
      is.maleYes = "Male",
      age.c = "Age",
      "commitmentYes:casualtiesYes" = "Commitment x Casualties",
      "commitmentYes:dyadYes" = "Commitment x Dyad"
    )
  ) +
  xlab("") + ylab("") + theme_bw() +
  scale_colour_grey(
    start = .2,
    end = .8,
    name = "Mediators",
    labels = c('Morality', 'Interest', 'Reputation', 'Identity')
  ) +
  scale_fill_grey(
    start = .2,
    end = .8,
    name = "Mediators",
    labels = c('Morality', 'Interest', 'Reputation', 'Identity')
  ) 

# Mediation Analysis
print("Mediation Model - Full Dataset")
# See Online Appendix F
process(data = data, 
        y = "defend",
        x = "commitment.b",
        m = c("morl","intr", "rept", "idnt"),
        cov = c("ideo.c", "news.c", "cybr.know.c", "polt.know.c", "coop.c", "age.c", "is.male.b"),
        effsize = 1,
        modelbt = 1,
        contrast = 1,
        seed = 1234,
        model = 4)

# See Online Appendix H

linr.morl = lm(morl ~ commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.morl.defd = lm(defend ~ morl + commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
morl.med = mediate(linr.morl, linr.morl.defd, treat = "commitment", mediator = "morl", robustSE = TRUE, sims = 5000, seed = 1234)
morl.sens = medsens(morl.med, rho.by = 0.1, effect.type = "indirect", sims = 5000)
print("Morality Sensitivity Analysis")
summary(morl.sens)
plot(morl.sens, sens.par = "rho")

linr.intr = lm(intr ~ commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.intr.defd = lm(defend ~ intr + commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
intr.med = mediate(linr.intr, linr.intr.defd, treat = "commitment", mediator = "intr", robustSE = TRUE, sims = 5000, seed = 1234)
intr.sens = medsens(intr.med, rho.by = 0.1, effect.type = "indirect", sims = 5000)
print("Interest Sensitivity Analysis")
summary(intr.sens)
plot(intr.sens, sens.par = "rho")

linr.rept = lm(rept ~ commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
linr.rept.defd = lm(defend ~ rept + commitment + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)
rept.med = mediate(linr.rept, linr.rept.defd, treat = "commitment", mediator = "rept", robustSE = TRUE, sims = 5000, seed = 1234)
rept.sens = medsens(rept.med, rho.by = 0.1, effect.type = "indirect", sims = 5000)
print("Reputation Sensitivity Analysis")
summary(rept.sens)
plot(rept.sens, sens.par = "rho")

# See Online Appendix G

print("Mediataion Model - Democracts")
process(data = subset(data, ideo < 4), 
        y = "defend",
        x = "commitment.b",
        m = c("intr","rept", "morl", "idnt"),
        cov = c("news.c", "cybr.know.c", "polt.know.c", "coop.c", "age.c", "is.male.b"),
        effsize = 1,
        modelbt = 1,
        contrast = 1,
        seed = 1234,
        model = 4)

print("Mediataion Model - Republican")
process(data = subset(data, ideo > 4), 
        y = "defend",
        x = "commitment.b",
        m = c("intr","rept", "morl", "idnt"),
        cov = c("news.c", "cybr.know.c", "polt.know.c", "coop.c", "age.c", "is.male.b"),
        effsize = 1,
        modelbt = 1,
        contrast = 1,
        seed = 1234,
        model = 4)

print("Mediataion Model - Independent")
process(data = subset(data, ideo == 4), 
        y = "defend",
        x = "commitment.b",
        m = c("intr","rept", "morl", "idnt"),
        cov = c("news.c", "cybr.know.c", "polt.know.c", "coop.c", "age.c", "is.male.b"),
        effsize = 1,
        modelbt = 1,
        contrast = 1,
        seed = 1234,
        model = 4)

# Robustness Checks

linr.cyber.treat = lm(resp.cybr.s ~ commitment + casualties + dyad, data = data)
linr.cyber.media = lm(resp.cybr.s ~ commitment + casualties + dyad + morl + intr + rept + idnt, data = data)
linr.cyber.full = lm(resp.cybr.s ~ commitment + casualties + dyad + morl + intr + rept + idnt + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)

# See Online Appendix C

print("Cyber - Treatment Model")
summary(linr.cyber.treat)
print("Cyber - Treatment + Mediator Model")
summary(linr.cyber.media)
print("Cyber - Full Model")
summary(linr.cyber.full)

linr.conv.treat = lm(resp.conv.s ~ commitment + casualties + dyad, data = data)
linr.conv.media = lm(resp.conv.s ~ commitment + casualties + dyad + morl + intr + rept + idnt, data = data)
linr.conv.full = lm(resp.conv.s ~ commitment + casualties + dyad + morl + intr + rept + idnt + ideo.c + news.c + cybr.know.c + polt.know.c + coop.c + is.male + age.c, data = data)

# See Online Appendix D

print("Conventional - Treatment Model")
summary(linr.conv.treat)
print("Conventional - Treatment + Mediator Model")
summary(linr.conv.media)
print("Conventional - Full Model")
summary(linr.conv.full)

# Balance Checks 

# See Online Appendix A
print("Commitment Balance Checks")
MatchBalance(as.numeric(commitment) - 1 ~ ideo + news + cybr.know + polt.know + coop + is.male + age, data = data, , nboots = 5000)
print("Casualties Checks")
MatchBalance(as.numeric(casualties) - 1 ~ ideo + news + cybr.know + polt.know + coop + is.male + age, data = data, , nboots = 5000)
print("Dyad Balance Checks")
MatchBalance(as.numeric(dyad) - 1 ~ ideo + news + cybr.know + polt.know + coop + is.male + age, data = data, nboots = 5000)
